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Abstract 

It is well known from numerous experiments that nuclear multifragmentation 
is a dominating mechanism for production of intermediate-mass fragments 
in nucleus-nucleus collisions at energies above 100 A MeV. In this paper we 
investigate the validity and performance of the Fermi break-up model and 
the statistical multifragmentation model implemented as parts of the Geant4 
toolkit. We study the impact of violent nuclear disintegration reactions on 
the depth-dose profiles and yields of secondary fragments for beams of light 
and medium-weight nuclei propagating in extended media. Implications for 
ion-beam cancer therapy and shielding from cosmic radiation are discussed. 
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1. Introduction 

The passage of energetic nuclei through various materials is a subject of 
experimental and theoretical investigations for several decades. In partic- 
ular, the energy deposition by light nuclei in living tissues is of a primary 
importance for ion-beam cancer therapy [H, [2j • Several Monte Carloparticle 
transport codes like SHIELD-HItJIl PHITS FLUKA 0,3, MC- 



NPX @, [h} and Geant4 [Til, E, HH are used to simulate the interaction 



* Corresponding author: pshenich@fias.uni-frankfurt.de 



Preprint submitted to Elsevier 



November 10, 2009 



of protons and nuclei with homogeneous phantoms made of tissue-like mate- 
rials, voxelized phantoms based on 3D CT images of patient's body, as well 
as with beam-line elements, e.g^ nozzles and collimators. 



As demonstrated in Refs. [15|, [16|, [17|, uM by calculations with the Monte 



Carlo Model for H eavy -Ion Therapy (MCHIT), the recent versions of the 



Geant4 toolkit [19|, [20| can be successfully used for simulations of carbon- 
ion therapy. However, until now the projectile nuclei heavier than carbon 
have received much less attention in studies with Geant4. Recently we have 
performed a comparative study of the depth-dose distributions for beams of 



various light and medium-weight nuclei in water using the Geant4 toolkit [21 



A proper description of physical and biological processes induced by light 
and heavy nuclei in extended media is important also for present and future 
exploration of space. This includes the d esig n of shielding elements to be 



used during future interplanetary missions [22], |23[ , reliable estimations of the 



doses acquired by astronauts in long space flights, evaluations of radiation 



conditions at Moon and other space colonies. As shown in Ref. [24J, the 
Galactic Cosmic Rays (GCR) present one of the main dangers for long-term 
human activities in space. The GCR present a very broad spectrum of nuclei, 
from hydrogen to iron, with kinetic energies from 10 MeV to 100 GeV per 
nucleon. The energy distribution has a maximum around 1 GeV per nucleon, 



however, the tails can reach 1 TeV per nucleon. It was demonstrated [24 
that despite of their small fraction in GCR (~ 10~ 3 ), the biological impact 
of heavy nuclei like Fe on living tissues is very strong, comparable with the 
impact of protons. Therefore, it is important to extend the capabilities of 
particle transport codes to the domain of heavy projectile nuclei, as it was 
done recently, for example, for the PHITS 0, and MCNPX ^ codes. 

During propagation of high-energy nuclei in extended media a certain 
fraction of them undergoes inelastic interactions with target nuclei. It is 
quite common to assume that this process consists of three stages: (1) the 
initial fast stage of reaction when nuclei interact strongly with each other; (2) 
preequilibrium stage when fast particles leave a highly excited nuclear system; 
and (3) de-excitation of the equilibrated nuclear residues by evaporation of 
nucleons and light clusters or by breaks up into several fragments. 

Most dynamical models used to describe the initial stage of nucleus- 
nucleus collisions come to the conclusion that after a time interval of a few 
tens of fm/c fast particles escape, and the evolution of the remaining nu- 
clear system changes its character. Because of intensive interaction between 
nucleons the residual nuclei evolve toward statistical equilibrium. If the ex- 
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citation energy is low, the hot compound nucleus undergoes de-excitation by 
evaporating nucleons and light clusters. At higher excitation energies the hot 
residual nucleus expands and breaks up into hot primary fragments, which 
later on undergo de-excitation by evaporating nucleons and light clusters. 
It is interesting to note that a wide variety of secondary nuclei, including 
exotic nuclei, e.g. hypernuclei, can be created in the Earth's atmosphere by 
energetic cosmic nuclei. However, it is very difficult to detect these products 
of nuclear reactions. 

In our studies with Geant4 15j, [ly, |17J the initial stage was simulated by_ 
means of the Light-Ion Binary cascade model (G4BinaryLightIonReaction) [27 
This model treats a nucleus-nucleus collision as a sequence of individual 
nucleon-nucleon collisions in the region where the colliding nuclei overlap 
with each other. The production of secondary hadrons, e.g. pions, in nucleon- 
nucleon collisions is also taken into account. Since this model is a straight- 



forward extension of the hadron- nucleus Binary cascade model 28J, it does 
not include collective phenomena such as compression of nuclear matter or 
collective flow of secondary products. Due to these simplifications the model 
developers suggested that it is relevant to nucleus-nucleus collisions where 
at least one of the colliding nuclei is not much heavier than carbon nucleus. 
In general, the Monte Carlo simulations for ion-therapy with beams of light 
nuclei 15|, ll_6|, [17J propagating in tissue-like materials are within the scope of 
the Binary cascade model, as mostly H, C, and O nuclei are involved. The 
Light-Ion Binary cascade model is recommended to be applied for projectile 
nuclei with energies between 50 A MeV and 10 A GeV [27|. The primary 
excitation energy of a nuclear residue left after the cascade stage of reaction 
is defined by the numbers of excitons,i.e. particles (captured nucleons) and 
holes (nucleons removed from the nuclear core). 

The Wilson abrasion model (G4WilsonAbrasionModel) is available as a 



part of Geant4 toolkit 29|, |30| as an alternative approach for simulating 



inelastic nucleus-nucleus collisions. It is dealing with calculating volumes of 
nuclear residues (spectators) and their excitation energies rather than with 
detailed consideration of nucleon-nucleon collisions in the overlap (fireball) 
zone. For the sake of simplicity the angular distribution of the nucleons 
abraded from a nucleus is assumed to be isotropic in the reference frame of 
this nucleus. The momentum distributions of such nucleons are simulated 
according to a phenomenological expression with parameters depending on 



the value of Fermi momentum of the corresponding nucleus |30J ■ The model 



has no restrictions on the masses of colliding nuclei, but the production of 
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secondary hadrons, e.g. mesons, in the overlap zone is ignored. Therefore, 
the predictions of the Wilson abrasion model become less and less accurate 
with increasing beam energy above the pion production threshold. This 
reduces the energy domain of applicability of the Wilson abrasion model as 
compared to the Light-Ion Binary cascade model. 

In the abrasion model the excitation energy of nuclear residues comes 
from two sources: (1) from excess of their surface energy associated with their 
non-spherical shape and (2) from the energy deposited by abraded nucleons 
passing through them. As well established, several de-excitation mechanisms 
are responsible for de-excitation of nuclear residues: emission of nucleons and 
light clusters, nuclear fission and multifragmentation. As follows from the 



statistical model analysis [311 ]. evaporation dominates at low excitation en- 
ergies, less than 3 MeV per nucleon. It is characterized by modest changes 
in the charge and mass of the initial residual nucleus, as only few neutrons, 
protons and alphas are emitted in a relatively long time scale. Nuclear fission 
is a dominant de-excitation mechanism only for very heavy nuclear residues 
like Th or U. The multifragment break-up of nuclei or multifragmentation 
becomes important at excitation energies above 3 MeV per nucleon and is 
characterized by a relatively short time scale of about 100 fm/c [311 ] . As es- 
tablished by many experiments, this is a dominating process for production 
of intermediate mass fragments. In recent years it has been demonstrated 
by numerous calculations that this process can be well described within the 
statistical approaches such as Fermi break-up model for light nuclei and Sta- 



tistical Multifragmentation Model (SMM) for medium and heavy nuclei 31 



In Ref. 32j, the contribution of multifragmentation to the total reaction cross 
section was estimated at the level of 10-15% for hadron-nucleus collisions, and 
about twice as large for nucleus-nucleus collisions. 

To the best of our knowledge, only the SHIELD-HIT code and the Geant4 
toolkit take into account multifragmentation of medium-weight and heavy 
nuclei in transport calculations for extended media. Due to radical changes 
in mass and charge of the projectile nucleus following multifragmentation, 
the ionization energy loss changes dramatically. Therefore, inclusion of nu- 
clear multifragmentation in calculations of ion energy deposition in extended 
media will lead to changes in depth-dose distributions as well as in yields of 
secondary fragments produced in nuclear fragmentation reactions. We be- 
lieve that the role of violent disintegration of highly excited nuclei should 
be carefully evaluated for a set of practical applications which deal with en- 
ergetic ion transport in extended media, especially for cancer therapy and 
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cosmic radiation protection. 

The paper is organized as follows. First, we give a brief description of the 
Fermi break-up model for decay of light nuclei (G4FermiBreakUp) in Sec. [21 
and, second, the statistical multifragmentation model (G4StatMF) in Sec. [5j 
Both models are parts of the Geant4 toolkit. Then we report the results 
of stand-alone tests of these models. Validation checks of G4FermiBreakUp 
and G4StatMF are made in Sec. Eland Sec. El respectively. We simulate mul- 
tifragment break-up of several hot nuclear systems with specified excitation 
energies. The results are compared with predictions of the FORTRAN-77 



version of the SMM [31[, which also includes the Fermi break-up model. In 
Sec. |H we calculate the depth-dose distributions for carbon ions in water, 
as well as the depth-yield distributions of secondary fragments produced by 
such nuclei. The calculational results are compared with corresponding ex- 
perimental data. The depth-dose distributions for iron nuclei in water are 
also calculated and compared with experimental data in Seed Finally, in 
the same section we demonstrate the impact of multifragmentation reactions 
on charge distributions of secondary fragments produced by iron nuclei in 
polyethylene. Our conclusions are formulated in Sec. [HI 

2. Fermi break-up model for violent fragmentation of light nuclei 

For a light nucleus with a mass number Aq and a charge Zq (in the fol- 
lowing we assume A < 16) even a relatively small excitation energy may 
be comparable to its binding energy. In this case we assume that the explo- 
sive decay of the excited nucleus into several smaller clusters is the principal 
mechanism of de-excitation. To describe this process we use a model which is 
similar to the famous Fermi break-up model for multiple particle production 



in proton-proton collisions [33(. Later on it was extended to the multifrag- 



ment break-up of highly excited nuclei, see, e.g., Refs. [34J, [35| and references 
therein. It is assumed that the excited nucleus decays simultaneously into 
cold or slightly excited fragments, which have lifetimes longer than the decay 
time, estimated as about 100 fm/c. The break-up configuration is character- 
ized by some freeze-out volume V, where the produced fragments are placed. 
The volume available for the translational motion, the so-called free volume 
Vf, is smaller than the freeze-out volume, at least, by the proper volume 
of the fragments, Vf = V — Vb, where Vo = A /p is the initial volume of 
the nucleus at normal nuclear density po = 0.15 fm -3 . Below the following 
parameterization is used: Vf = kVq, while k is a model parameter of order 1. 
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The masses of fragments in their ground and lowest excited states were taken 



from nuclear data tables 36 



We consider all possible break-up channels, which satisfy the mass num- 
ber, charge, energy and momenta conservations, and take into account the 
competition between these channels. The central assumption of the model 
is that the probabilities of different break-up channels are proportional to 
their microcanonical weights. In other words the probability of an individ- 
ual break-up channel containing n particles with masses mi (z = 1, • • • , n) is 
proportional to its phase space volume [33j, |34|, |35|, l37j |: 



w: 
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where M = YH=i m i ^ s the total mass of the fragments, S = YYi=i(^ s i + -0 
is the spin degeneracy factor (sj is the i-th particle spin), G = rij=i n j' * s 
the particle identity factor (rij is the number of particles of kind j). Ekin is 
the total kinetic energy of particles at infinity which is related to the nucleus 
excitation energy E\ z as 



E ki n = E* AoZo + M c 2 - ^m iC 2 . 



(2) 



i=l 



M is mass of the decaying nucleus, is the Coulomb interaction energy 



between fragments given in the Wigner-Seitz approximation [31 
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where Ai, Zi are mass numbers and charges of produced particles. In our 
calculations we have included fragments in all excited states, which are stable 
with respect to the nucleon emission, as well as long-lived unstable nuclei 5 He, 
5 Li, 8 Be, 9 B, which decay at the later stage of the reaction 35]. 

The number of channels included in our calculations was about 10 3 for 
the 16 O nucleus and ~ 2 ■ 10 2 for the 12 C. The Coulomb expansion stage 
was not considered explicitly for such light systems. The momentum distri- 
butions of final products were obtained by the random generation over the 
whole accessible phase space, determined by the total kinetic energy, Eq. (j2J), 
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taking into account the energy and momentum conservation. For the calcu- 
lations of the available phase-space a very effective algorithm proposed by 
G.I. Kopylov 



38, 39 was used. 



The Fermi break-up model in combination with the intranuclear cascade 
model well describes various experimental data on disintegration of light nu- 



clei by energetic protons with kinetic energies above 100 MeV |3ll . l35l . |4Q| . |41 
However, as pointed out in Ref. 



35 



the model becomes less accurate below 
100 MeV when the emission channels with only few fragments dominate, in 
particular, in the case of 12 C and 16 targets. The description of such two- 
and three-body decay channels can be improved by considering alternative 
reaction mechanisms, e.g. the direct knock-out of light clusters, which do 
not requires the formation of a compound nucleus. 



3. Validation of Fermi break-up model of Geant4 

The Fermi break-model for decay of light nuclei was implemented as a 
FORTRAN- 77 code, see, e.g., Ref. | 3lU35lj . Later the same model was imple- 
mented in C++ by Vicente Lara [42, |43| and became a part of the Geant4 
toolkit as the G4FermiBreakUp class. The G4FermiBreakUp has a method 
BreakltUp which is applicable to Geant4 objects of the type G4Fragment 
representing excited nuclei. The BreakltUp results in a G4Fragment Vector 
which contains a set of G4Fragment objects. We can investigate now whether 
the results of the both implementations of the Fermi break-up model give con- 
sistent results. Before running tests we have corrected the table of energy 
levels of excited light nuclei used by G4FermiBreakUp. 

3.1. Average multiplicities of nuclear fragments calculated with the Fermi 
break-up model 

The average multiplicity of nuclear fragments created in decays of ex- 
cited nuclei is an important characteristic of the fragmentation process which 
shows clearly its violent nature. The average multiplicities calculated with 
the FORTRAN and C++ implementations of the Fermi break-up model are 
plotted in Fig. [T] as functions of excitation energy for decays of 12 C, 13 C, 12 N 
and 13 N. Such highly-excited nuclei are abundantly produced, in particular, 
in interactions of therapeutic carbon-ion beams with tissue-like materials. 
These nuclei represent either excited projectile carbon nuclei or fragments of 
oxygen nuclei from the media. As seen from the figure, the predictions of 
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two codes for average fragment multiplicities in decays of such nuclei are in 
very good agreement. 




Excitation energy (MeV/nucleon) Excitation energy (MeV/nucleon) 



Figure 1: Color online. Average multiplicities of nuclear fragments created in decays of 
excited 12 C, 13 C, 12 N and 13 N nuclei as functions of their excitation energies. The results 
of the G4FermiBreakUp of Geant4 are shown by histograms, the results of the Fermi 
break-up model of the SMM code are shown by points. 

3.2. Average multiplicities of hydrogen, helium, lithium and beryllium frag- 
ments 

The average multiplicities of hydrogen, helium, lithium and beryllium 
fragments created in decays of 12 C are shown in Fig. [2] as functions of 12 C 
excitation energy. One can see that the H and He isotopes are most abundant 
decay products at all excitation energies. As seen in Figs. [1] and [21 both 
implementations of the Fermi break-up model predict a peak at excitation 
energies of 1-2 MeV/nucleon corresponding to the decay of 12 C into three 



S 



a-particles. At excitations above 5-6 MeV/nucleon the nuclei decay mostly 
into protons and neutrons, but on average there is also one a-particle among 
decay products. The production of Li and Be fragments reaches maximum at 
excitation energies of 7-8 and 5-6 MeV/nucleon, respectively. The differences 
between decays of 12 C, 13 C, 12 N and 13 N nuclei reveal themselves only at low 
excitations, e.g. below 2 MeV/nucleon. At 10 MeV/nucleon all nuclei decay 
into 4-5 fragments, on average, as predicted by the both codes. 




Excitation energy (MeV/nucleon) Excitation energy (MeV/nucleon) 

Figure 2: Color online. Average multiplicities of H, He, Li and Be fragments created in 
decays of excited 12 C as functions of its excitation energy. The notations are the same as 
in Fig. [TJ 

In general, the results obtained with two implementations of the Fermi 
break-up model agree well at low and high excitations for all fragments, from 
H to Be. Some minor discrepancies in results are seen only for the yields of He 
and Be nuclei between 2.5 and 6 MeV/nucleon. We attribute these deviations 
to different nuclear mass tables used by the considered codes. 
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4. Production of secondary fragments by light nuclei in tissue-like 
media 



The depth-profiles of average linear energy deposition calculated per beam 
particle for 200 and 400 A MeV 12 C in water are shown in Figs. [3] and HI 
They were obtained with the Light-Ion Binary cascade model and the Wilson 
abrasion model, respectively. The calculations were performed for a water 
phantom of 40 x 15 x 15 cm, which was divided into 0.25 mm slabs. The 
energy deposition was calculated for each of the slabs. At the end of the 
run it was divided by the slab thickness and by the number of projectile 
nuclei to estimate the average linear energy transfer per beam particle. A 
Gaussian beam cross section of 5 mm FWHM was assumed. The beam 
energy distribution was taken as a Gaussian with the FWHM of 0.2% of the 
mean energy. 

Experimental data of Ref. 44J for relative ionization in water for 200 and 
400 A MeV 12 C beams were rescaled to the calculated absolute dose values 
at zero depth. These data are plotted in Figs. [3J and 0] and compared to the 
calculations using both nucleus-nucleus collision models mentioned above. In 
both cases the MCHIT model successfully describes the positions of the Bragg 
peak and the peak-to-plateau ratios of the depth-dose curves for 200 and 
400 A MeV 12 C ions in water. We used the set of standard electromagnetic 
models of Geant4 for calculating ionization energy losses of carbon nuclei and 
their fragments. For better description of the positions of the Bragg peaks 
the average ionization potential for water molecule was set to 81 eV, which 
is in the range of uncertainty quoted for this quantity 45 . 

The normalized yields of primary carbon nuclei and secondary nuclear 
fragments (hydrogen, helium, lithium, beryllium and boron nuclei) were also 
calculated with MCHIT and plotted in Figs. [3J and HI In thes e fig ures the 
calculated yields are compared with the measurements of Ref. [44j. As the 
yields were measured by detecting fragments within a 10° cone surrounding 



the beam axis, as explained in Ref. [44| , corresponding acceptance cuts were 
introduced in calculations. 

Both the Light-Ion Binary cascade and Wilson abrasion models describe 
rather well the observed attenuation of the carbon beam in water at 200 and 
400 A MeV. The measured and calculated carbon yields correspond to the 
mixture of primary 12 C and secondary carbon nuclei like 9 C, 10 C and n C. 
The last two nuclides ( 10 C and 11 C) are relevant to PET monitoring of dose 



distribution in carbon-ion therapy. In Ref. [16[ their yields were calculated 
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Figure 3: Depth-dose distributions and normalized yields of secondary fragments for 200 
and 400 A MeV 12 C in water calculated with MCHIT involving the Light-Ion Binary 
cascade model: (a) the average linear energy deposition per beam particle; the yields of 
nuclear fragments per beam particle (b) for carbon, helium and hydrogen fragments; (c) 
for boron, beryllium and lithium fragments. The yields of secondary fragments calculated 
with G4FermiBrcakUp applied after the cascade model are shown by thick lines. The 
same yields, but calculated with the evaporation model applied after the Light-Ion Binary 
cascade model are shown by thin lines. Experimental data from Ref. [3] are shown by 
various symbols. 
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Figure 4: Dcpth-dosc distributions and normalized yields of secondary fragments for 200 
and 400 A MeV 12 C in water calculated with MCHIT involving the Wilson abrasion model. 
The notations are the same as in Fig. [3] 
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with the MCHIT model and compared with available experimental data. As 
follows from the data and calculations, more than 40% of primary 200 A MeV 
12 C nuclei undergo fragmentation on their way to the Bragg peak, while at 
400 A MeV this fraction exceeds 70%. 

The production of H, He, Li, Be and B nuclei is characterized by gradually 
increasing yields until primary ions stop at the depth, where the Bragg peak 
is located. Secondary fragments propagate further beyond the distal edge of 
the Bragg peak, and their yields decrease due to secondary fragmentation 
reactions on nuclei in the media. The boron yield drops fastest since B has 
the largest charge and inelastic cross section among the elements produced 
by 12 C ions. On the contrary, H and He fragments are less attenuated and 
most of them leave the water phantom of 40 cm depth. 

By inspecting Figs. [3] and H] one can compare the Light-Ion Binary cascade 
model and the Wilson abrasion model with respect to their ability to describe 
secondary fragments. The yields of secondary fragments calculated with the 
evaporation model instead of the Fermi break-up model are also shown for 
comparison. 

The cascade model underestimates the yields of hydrogen fragments, 
while the abrasion model successfully describes them when the Fermi break- 
up model is used instead of the evaporation model. Both nucleus-nucleus col- 
lision models underestimate He yields, and the inclusion of G4FermiBreakUp 
does not improve this situation. This deficiency is apparently caused by ne- 
glecting the alpha-clustered structure of 12 C, which may enhance the pro- 
duction of a-particles (see also Ref. [401]). 

As seen from Fig. [3], the cascade model reasonably describes the B yields, 
but fails to reproduce the production of Li and Be. When the evaporation 
model is used instead of G4FermiBreakUp, the agreement with Be data is 
achieved, but the production of Li is strongly underestimated. The abrasion 
model gives more accurate predictions for the yields of Li and Be fragments, 
but only in combination with the Fermi break-up model. When instead the 
evaporation model is used, the Li and Be yields are underestimated by large 
factors. At the same time both de-excitation models underestimate the B 
yields. This can be explained either by the failure of the abrasion model to 
describe peripheral nucleus-nucleus collisions, leading to low-excited residual 
nuclei, or by the inaccuracy of the Fermi break-up model in describing the 
two-body decays of nuclei with lowest excitations. 

In summary, our analysis shows that no one of the considered models 
describes the production of secondary fragments perfectly. Generally, the 
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fragment yields are somewhat better described by the abrasion model as 
compared with the cascade model. However, the former model still under- 
estimates the yields of He and B fragments. The inclusion of the Fermi 
break-up model at the de-excitation stage of nuclear reaction improves the 
description of H, Li and Be yields, but does not help to remove discrepancies 
for B fragments. This is an important issue for Monte Carlo calculations for 



carbon-ion therapy. Indeed, as shown recently [46(, the dose and radiation 
quality at the center of a broad carbon-ion beam formed by using a passive 
beam delivery system are influenced by such secondary fragments from the 
off-center region. 

5. Break-up of medium- weight and heavy nuclei simulated with 
SMM model 

The Statistical Multifragmentation Model (SMM) is based on the assump- 
tion of statistical equilibrium between produced fragments in a low-density 



freeze-out volume [31[. We believe that the chemical composition, i.e. masses 
and charges of primary fragments, are fixed at this stage. However, the 
fragments can still interact with other nuclear species via the Coulomb and 
nuclear mean fields. Hence their energies and densities may be affected by 
these residual interactions. 

An advantage of the model is that all breakup channels to nucleons and 
excited fragments are considered within the same statistical framework, and, 
in particular, the formation of a compound nucleus is included as one of the 
channels. This allows for a smooth transition from the decay via evaporation 



and fission at low excitation energies [47| to the multifragmentation at high 
excitations. 



In the microcanonical treatment |3ll . l48j the mass, charge, momentum 



and energy of the system are strictly fixed. It is also assumed that the 
primary fragments are formed in the expanded volume V > Vq, where Vq is 
the volume at normal nuclear density po — 0.15 fm -3 . 

In accordance with the statistical hypothesis the probability of the decay 
channel j is given by the statistical weight Wj oc expS^, where Sj is the 
entropy of the system in channel j which is a function of the excitation 
energy E x , mass number A , charge Z and other global parameters of the 
source. After formation in the freeze-out volume, the fragments propagate 
independently in their mutual Coulomb field and undergo secondary decays. 
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De-excitation of the hot primary fragments proceeds via evaporation, fission, 



or Fermi-breakup [37 



In the SMM light fragments with mass number A < 4 and charge Z < 
2 are considered as structureless particles (nuclear gas) with their masses 
and spins taken from the nuclear data tables. Only translational degrees of 
freedom of these particles contribute to the entropy of the system. Fragments 
with A > 4 are treated as heated drops of nuclear liquid, and their individual 
free energies Faz are parameterized according to the liquid drop model: 



Flz + Flz + E C AZ + E S A T, (4) 



where terms in the r.h.s correspond to the bulk, surface, Coulomb and sym- 
metry energy. In this expression F AZ = (—Wo — T 2 /e )A is the bulk en- 
ergy term including the contribution of internal excitations controlled by the 
level-density parameter eo, and Wq = 16 MeV is the binding energy of infinite 
nuclear matter. F% z = B A 2/3 ((T 2 - T 2 )/(T 6 2 + T 2 )) 5/4 is the surface energy 
term, where Bq = 18 MeV is the surface coefficient at T = 0, and T c = 
18 MeV is the critical temperature of infinite nuclear matter. The Coulomb 
energy of individual fragments is calculated as E AZ = f ^1/3 c (p) > where e 
is the proton charge, r =1.17 fm, and the last factor, c(p) = 1 — (p p / Ppo) 1 ^ 3 , 
where p p o is the normal proton density of nuclei p p o ~ j^Po- This factor 
describes the screening effect due to the presence of other fragments in the 
Wigner-Seitz approximation. The last term in Eq. (j4J) E s Az m = r )(A — 2Z) 2 /A 
is the symmetry energy term, where 7 = 25 MeV is the symmetry energy 
coefficient. These parameters are taken from the Bethe-Weizsacker formula 
and correspond to the isolated cold fragments with normal nuclear density. 
This assumption has been proven to be quite successful in many applications. 
However, a more realistic treatment of primary fragments in the freeze-out 
volume may require certain modifications of the liquid-drop parameters as 



suggested by experimental data [49|, |50|, |51 |. 



In the grand canonical (macrocanonical) version of the SMM [52J, after 
integrating out translational degrees of freedom, one can write the mean 
multiplicity of nuclear fragments with mass A and charge Z as 



(N A z) = g A z^A 3/2 ex P 



T 



±(F A z(T,p)-pA-vZ) 



Here g A z is the ground-state degeneracy factor of species (A,Z), At = 

1 /2 

(2-Kh 2 /ttinT) is the nucleon thermal wavelength, and pa 939 MeV is 
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the nucleon mass. Here Vf ~ kVq is the free volume available for the transla- 
tional motion of fragments and k is the model parameter which in principle 
can depend on the fragment multiplicity in the freeze-out volume 31|. The 
chemical potentials /i and v are found from the mass and charge constraints: 



(A,Z) 



AZ 



)A = A Q , 



(A,Z) 



AZ 



(5) 



Numerous comparisons of the SMM calculations with experimental data 
on thin targets show that gen erally the model describes data very well (see, 



e.g., Refs. [54j, |55|, |56|, |57|, |58|, |59|, [60] ). This demonstrates that the statis- 



tical approach with liquid-drop description of individual fragments provides 
adequate treatment of the multifragmentation process. 



6. Validation of statistical multifragmentation model of Geant4 



6.1. General remarks 

The statistical approach to multifragment break-up of hot nuclear systems 



outlined above was first formulated in Refs. [52|, |53j, and it is well-known now 
as a Copenhagen-Moscow model (for a review see Ref. 3JJ). For numerical 



calculations it was initially implemented as a FORTRAN-77 SMM code, 



sec 



e.g., Ref. (37 



Lara 43, 61 



Later the SMM was implemented in C++ by Vicente 
and became a part of the Geant4 toolkit as the G4StatMF 
class. G4StatMF has a method BreakltUp which is applicable to a Geant4 
object of type G4Fragment representing a nucleus in its ground or excited 
state. After application of this method a G4Fragment Vector is produced as 
an output, which consists of a set of G4Fragment objects. 

It is worthwhile to mention that the C++ implementation of the SMM 



has been developed following the FORTRAN version of the model [37j. After 



inspecting the source code of G4StatMF we came to the conclusion that the 



physical parameters of the model are basically the same as in Refs. [3ll . |48 
However, the numerical methods involved to find, e.g. the temperature and 
chemical potentials for a macrocanonical ensemble of nuclear fragments, are 
different in the FORTRAN and C++ implementations of the SMM. The 
FORTRAN SMM works with single precision floating point numbers, while 
the G4StatMF uses double precision floats. In view of these differences we 
performed a systematic comparison of numerical results delivered by the two 
codes. 
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The G4StatMF from the Geant4 toolkit of version 9.1 was used to simu- 
late the multifragment break-up of hot nuclear systems with mass and charge 
of 112 Sn and 208 Pb. In stand-alone tests 10 5 decay events were generated for 
each nuclear system and excitation energy of 3, 4, 5 and 8 MeV/nucleon. In 
these calculations the parameter k which defines the free volume for transla- 
tional motion of fragments was taken as a function of each event's fragment 
multiplicity as k(Mj), where Mf = ^ ^ -/V^. During the simulation several 
kinds of distributions of produced fragments were scored. Such distributions 
were compared with the results obtained with the FORTRAN SMM using 
the same calculational parameters. 

6.2. Fragment mass distributions 

At the beginning of our study severe discrepancies between the codes were 
found for the mass distributions of produced fragments. The disagreements 
were already seen for the decay of 112 Sn, but they became more pronounced 
when we considered nuclear systems far from the stability line, e.g. 55 Sn and 
140 Sn. Such exotic proton- and neutron-rich systems cannot be produced 
in collisions of stable nuclei, but we consider them only for the verification 
of the code. Our simulations of multifragment break-up of such exotic nu- 
clear systems lead us to the conclusion that the calculation of the fragments' 
symmetry energy in G4StatMF of version 9.1 was wrong. The error was 
localized and corrected. After that the results of the two implementations of 
the multifragmentation model turned out to be very close. For example, the 
ensemble's average temperatures in decays of 112 Sn and 140 Sn differed by less 
than 0.2 MeV instead of 1 MeV before the corrections. In addition, several 
changes were introduced to improve the stability of G4StatMF when applied 
to proton- and neutron-rich systems like 23 0, 55 Sn and 50 Ne, 140 Sn. After 
these improvements to the G4StatMF were introduced, the corresponding 
distributions were calculated again in order to see the effect of updates. 

The results for 112 Sn are shown in Fig. [5]for excitation energies of 3 and 4 
MeV/nucleon where competition of explosive multifragment break-up and 
evaporation from the compound nucleus is taking place. This is clearly 
reflected in the mass distributions by the presence of heavy evaporation 
residues with A ~ 110. However, as clearly seen at 3 and 4 MeV/nucleon, the 
G4StatMF predicts much stronger compound nucleus peak than the standard 
SMM. Apparently this can be explained by a larger fraction of evaporation 
events produced by G4StatMF. At higher excitation energies the mass distri- 
butions produced with both codes agree quite well. At 5 and 8 MeV/nucleon 
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Figure 5: Color online. Mass distributions of nuclear fragments created after the decay 
of excited 112 Sn with excitation energies of 3, 4, 5 and 8 MeV/nucleon. The results of 
the modified G4StatMF of Geant4 are given by histograms, the results of the FORTRAN 
SMM - by open circles. 



the peaks associated with evaporation residues completely disappear and 
multifragment break-up becomes a dominating channel. It is seen also, that 
the production of intermediate mass fragments (IMF) with 4 < A < 20 is 
enhanced. At such high excitations the results of both codes are in perfect 
agreement both in shape and in absolute values. 

In Fig. [6] the results for a heavy system, 208 Pb, are shown for excitation 
energies of 3, 4, 5 and 8 MeV/nucleon. At 3 MeV/nucleon the contribution of 
fission-like events is clearly visible. Such events are characterized by creation 
of two fragments with comparable mass numbers of A ~ 100. Both SMM 
implementations predict such a behavior. At higher excitations the peak 
from symmetric fission-like events disappears. As seen from Fig. [6j the mass 
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distributions of fragments calculated with G4StatMF and FORTRAN SMM 
are in good agreement both in shape and in absolute values. 




Figure 6: Color online. Same as in Fig. [5] but for the decay of excited 208 Pb nucleus. 

As additional test for the consistency of the two SMM implementations 
one can calculate the average charge of the fragment as a function of the 
fragment mass, (Z)(A). It is expected that for each fragment mass the ratio 
(Z)/A is close to the ratio for the initial hot system. We found that the 
results of G4StatMF and the FORTRAN version of SMM follow this trend. 
Since the results of the two codes completely coincide, we do not present the 
corresponding plots here. 

6.3. Fragment charge distributions 

Further validation of G4StatMF is possible with experimental data on 
low-energy (several tens of MeV/nucleon) nucleus-nucleus collisions. This 
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energy regime allows colliding nuclei to overcome the Coulomb barrier and 
create a single combined system. There exist dedicated experimental tech- 
niques 62|,|63j, which make possible to select central collision events, resulting 
in a single emitting source, and allow to estimate the source's charge, mass 
and excitation energy. 

In particular, differential charge multiplicity distributions were measured 
in Ref. 62| for central 129 Xe+ nat Sn collisions at 32 A MeV. The event selection 



procedure adopted in Ref. [62[ required detection of a significant fraction 
(80%) of the total charge and of the initial linear momentum of colliding 
nuclei. Therefore, one can estimate the mass and charge of the emitting 
source as A = 198 and Z = 83 roughly corresponding to a very hot 198 Bi 
nucleus. Events containing less than three fragments with charge Z > 5 
were rejected to favor events resulting from the multifragmentation process. 
Differential charge multiplicity distributions for decay of 198 Bi were calculated 
with G4StatMF and plotted in Fig. [7J We have found that the experimental 
data are best described under the assumption that the 198 Bi system has the 
excitation energy of 5.9 MeV/nucleon. This value is lower compared to the 
available excitation energy per nucleon for a hypothetic compound system 



formed via a complete fusion of Xe and Sn nuclei [62J. This difference can 
be explained by the fact that the most energetic nucleons leave the system 
before the statistical equilibrium is reached, so that the average excitation 
energy per nucleon drops down. 

The fragment charge distributions were filled according to the event se- 



lection procedures adopted in the experiment [62J. In particular, following 
simulation of each individual event, the multiplicity of fragments with Z > 5 
in the event was calculated. Then the numbers of fragments of certain charge 
in the event were calculated and divided by this multiplicity value. After that 
the corresponding histogram shown in Fig. [7J was filled with these numbers. 

Distribution of hot fragments produced in multifragmentation events are 
shown by the dotted histogram in Fig. [7J As one can see, they significantly 
differ from the experimental data. We have found that the agreement with 
data is improved significantly if we account for secondary de-excitation of 
these fragments. This is especially important for heavy fragments with Z > 
30. The final distributions calculated with G4StaMF agree well with the 



experimental data of Ref. [621 ]. as well as with the distributions calculated 
with the FORTRAN SMM. 



20 




Figure 7: Color online. Fragment charge distributions (normalized to each event's Z > 5 
multiplicity) in decay of hot 198 Bi-like system representing the equilibrated source formed 
in central 129 Xc+ nat Sn collisions. The results of modified G4StatMF with and without 
accounting for secondary de-excitations of produced fragments are given by solid- and 
dashed-line histograms, respectively. The results of the FORTRAN SMM arc presented 
by open circles. Triangles represent experimental data [62]. 



6.4- IMF multiplicity distributions 

Now one can calculate the multiplicity distributions of intermediate mass 
fragments (IMF, 4 < A < 20) produced at the break-up of a 198 Bi-like nucleus 
created in central 129 Xe+ nat Sn collisions. Taking into account secondary 
decays the average IMF multiplicity calculated with the G4StatMF is 4.8. 
This value is slightly higher than the experimental value of 4.3 reported in 



Ref. [62J. The IMF multiplicity distribution calculated with G4StatMF is 



shown in Fig. [S] together with the distribution measured in Ref. [62[ . 

As seen from Fig. [HJ the general shape of the experimental distribution is 
well reproduced by G4StatMF, but the model distribution is shifted towards 
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Bi 5.9 MeV/nucI 
— G4StatMF 



A Data by Rivet et al. 




Figure 8: Color online. Multiplicity distributions of intermediate mass (Z > 5) fragments 
in decay of hot 198 Bi-like system representing the equilibrated source formed in central 
129 Xe+ nat Sn collisions. The results of modified G4StatMF are given by histograms. Tri- 
angles represent experimental data (6^ |. 



higher multiplicities. We attribute this deficiency to the uncertainties in the 
choice of the system size, which is smaller than the 198 Bi-source used in the 



simulations. As shown in Ref. [62j, the IMF multiplicity distributions are 
rather sensitive to the size of the emitting source. 

At the same time, the uncertainties in system size do not affect very much 
the fragment charge distributions shown in Fig. 0, where the experimental 
data are described well even without the exact specification of the system 
size in calculations. Indeed, as pointed out in Ref. the fragment charge 
distributions for hot nuclear systems are insensitive to the system size but 
are mostly governed by the system's excitation energy per nucleon. 
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7. Role of multifragmental ion in transport calculations of medium- 
weight nuclei in extended media 



After the validity of the G4StatMF has been confirmed, we can estimate 
the influence of multifragment break-up on ion transport in tissue-like mate- 
rials. It is evident that the yields of secondary fragments, especially of IMFs 
4 < A < 20, produced by beams of heavy nuclei in thick targets depend on 
the model used to simulate de-excitation of residual nuclei. Before consider- 
ing fragment yields, we examine more rough characteristic of ion transport 
through tissue-like materials, namely, the depth-dose distributions in water. 

7.1. Depth-dose distributions of 56 Fe beams in water 



The MCHIT model [18j described above has been extensively used to 



study ion propagation in tissue-like materials. Several examples of the depth 



dose distributions relevant to ion-beam cancer therapy can be found in Refs. [15, 



16|, 113) uM ■ Here we present results for heavier incident nuclei, Fe, which are 
relevant to the evaluation of cosmic radiation effects. 

The depth-dose distribution for 969.8 A MeV and 1087 A MeV 56 Fe nu- 
clei in water calculated with MCHIT are shown in Fig. [9] together with the 



experimental data of Refs. [26j and [64j, respectively. A water phantom of 



60 x 10 x 10 cm was divided into 0.5 mm slabs, and the average linear en- 
ergy deposition was calculated in each of the slabs. The calculations were 
performed for a Gaussian beam profile of 4 mm FWHM in the transverse 
plane to the beam direction. The energy spread of the beam was assumed 
to be Gaussian with the FWHM of 0.2% of the reported mean beam energy. 
Similar to calculations presented in Sec. HJ two models describing the ini- 
tial stage of inelastic nucleus-nucleus collisions, the Light-Ion Binary cascade 



model [27( and the Wilson abrasion model [29l . |30| were used in independent 
runs. 

De-excitation of hot nuclear residues formed in collisions of Fe ions with 
nuclei of the medium after emission of fast nucleons was treated by us- 
ing the BreakltUp method of the class G4ExcitationHandler. Within the 
G4ExcitationHandler several de-excitation models can be applied, namely, 
nucleon evaporation, photon emission, Fermi break-up and nuclear multi- 
fragmentation. The G4ExcitationHandler has public methods which set the 
domains of applicability of the de-excitation models in terms of the size and 
excitation energy of nuclear residues. However, since the Light-Ion Binary 
cascade model creates its own private instance of G4ExcitationHandler, such 
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Figure 9: Color online. Depth-dose distributions of (a) 969.8 A MeV and (b) 1087 A MeV 
56 Fe nuclei in water calculated with the MCHIT. The results of the Light-Ion Binary cas- 
cade model of Geant4 with and without multifragmcntation are shown by the dash-dotted 
and dotted histograms, respectively. The results of the Wilson abrasion model with and 
without multifragmcntation are shown by the solid and dashed histograms, respectively. 
The experimental data from Refs. 2(| 64 1 are shown by various points. 
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methods are currently not available to the user when this model is involved. 
By default, the multifragmentation is switched off, and we have changed 
the source code of the G4ExcitationHandler to activate it. On the contrary, 
The G4WilsonAbrasionModel has public methods which control the behav- 
ior of the G4ExcitationHandler. The distributions presented in Fig. [9] were 
calculated both with and without accounting for multifragment break-up at 
excitation energies > 3 MeV/nucleon. 

The probability of nuclear interactions for high-energy Fe nuclei in water 
is rather high. One can estimate the mean free path of the beam with 
respect to nuclear interactions, A nuc = l/(o" nuc n water (o)) ~ 17 cm, where 
0nuc ~ 2 barn is the total nuclear reaction cross section in Fe+O collisions 
and n watcr ( ) = 3 x 10 22 cm -3 is the density of oxygen nuclei in water. 

Nuclear interactions are reflected in the shape of the depth dose-distribution 
presented in Fig. [9j A globally descending trend is explained by the destruc- 
tion of beam particles in fragmentation reactions on nuclei of the medium. 
These reactions lead not only to reduction of the Bragg peak, but also to 
build-up of a tail beyond the peak due to deep penetration of light reaction 
products like protons or alphas. These trends have been alread y s een in 



MCHIT calculations for light and medium- weight nuclei in water [2l|. The 



non-trivial shape of the depth-dose curves shown in Fig. [9] is specific for 
deep-penetrating high-energy nuclear beams. 

Moreover, the ionization energy loss, proportional to the square of the 
nuclear charge drops dramatically due to nuclear fragmentation since {Z\ + 
Z 2 + ... + Zn) 2 > Z\ + Z\ + ...Zff. Nuclear evaporation usually leaves a 
heavy residue, while nuclear multifragmentation creates several fragments 
of comparable mass. In the latter case the total ionization energy loss of 
secondary fragments is reduced considerably compared to typical energy loss 
of evaporation residues. 

These considerations explain the trends seen in Fig. [9j in particular, why 
the linear energy deposition calculated with accounting for multifragmenta- 
tion is reduced compared to the calculations, where light particle evaporation 
is the only de-excitation process. However, as seen in Fig. [HI neglecting multi- 
fragment decays gives only very small changes in the depth-dose distributions 
calculated both with the cascade and abrasion models. This reflects the fact 
that the contribution of multifragmentation to the total reaction cross section 
is small. 

It should be stressed, that the calculations with the cascade model over- 
estimate the linear energy deposition from 56 Fe nuclei in water, while the 
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calculations with the abrasion model agree very well with the experimental 
data at two beam energies. We attribute this deficiency to the fact that 
the Light-Ion Binary cascade model underestimates either the fragmentation 
probability of heavy-nuclei, or the excitation energy of nuclear spectators 
produced after the fast stage of the reaction. 

However, as also seen from Fig. the results obtained with the abrasion 
model for 969.8 A MeV Fe underestimate the dose after the distal edge of the 
Bragg peak. This analysis demonstrates the necessity to improve physical 
models for the description of the initial stage of high-energy nucleus-nucleus 
collisions especially for medium-weight and heavy projectiles. 

7.2. Yields of secondary fragments produced by medium-weight nuclei in ex- 
tended media 

Further validation of the nucleus-nucleus collision models of the Geant4 
toolkit can be performed with the data on the yields of secondary fragments 
produced by nuclear beams in thin and thick targets. In particular, such 
data were presented in Refs. 0, [64| for 56 Fe nuclei transversing polyethylene 
(PE) targets. Calculations with the MCHIT model were performed for such 
polyethylene targets with areal densities of 1.94, 4.2 and 17 g/cm 2 . The 
yields of charged nuclear fragments leaving such targets at any angle were 
scored according to their charge. The results were divided by the number of 
primary beam nuclei to obtain normalized fragment yields which were plotted 



in Fig. [TU] together with the data of Refs. |5j, |64|. Again, the cascade and 
abrasion models were used in the MCHIT calculations. 

In Fig. [10] calculations with and without multifragment decays are pre- 
sented for 4.2 and 17 g/cm 2 polyethylene targets making possible to estimate 
the impact of explosive fragmentation on fragment yields predicted by the 
two nucleus-nucleus collision models. As seen from Figs. [TDT b) andfTOTc). the 
yields of intermediate mass fragments with Z = 6 — 9 can not be properly 
described by either of two models without taking into account multifragment 
decays. 

The general shape of Z-distributions is better reproduced by the cascade 
model in combination with G4StatMF. The abrasion model with or without 
multifragmentation underestimates the fragment yields in a broad range of 
charges from 10 to 22. Since the largest fragments contribute the main 
part to the total dose, the models should especially accurately estimate the 
production of fragments close in mass to the beam nuclei. As seen from 
Fig. fTOT b). both models describe the yields of Z = 25 and Z = 26 nuclei 
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Figure 10: Color online. Yields of charged nuclear fragments produced by 1050 A MeV 
and 998.8 A MeV 56 Fc nuclei in polyethylene targets of (a) 1.94, (b) 4.2 and (c) 17 g/cm 2 
areal densities. The MCHIT results obtained with the Light-Ion Binary cascade model of 
Geant4 with and without multifragmentation are shown by the dotted and dash-dotted 
histograms, with the Wilson abrasion model with and without multifragmentation are 
shown by the solid and dashed histograms, respectively. The experimental data from 
Refs. p, EH are shown by various points. 
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on 4.2 g/cm 2 target. However, the cascade model significantly overestimates 
such yields on a more thick target of 17 g/cm 2 , see Fig. ITOlf c). while the 
abrasion model is closer to experimental data for this case. This may explain 
why the abrasion model better describes the depth-dose distributions of Fe 
nuclei in thick targets made of light materials, as it was already demonstrated 
in Sec. 17.11 for water. 



8. Conclusions 

We have performed validation checks of the C++ versions of the Fermi 
break-up model (G4FermiBreakUp) and the Statistical Multifragmentation 
model (G4StatMF) of the Geant4 toolkit (version 9.1) by comparing their 
results for multifragmentation of light, medium-weight and heavy nuclei with 
results of the FORTRAN versions of the corresponding models and with 
available experimental data. Two versions of the Fermi break-up model agree 
well after the corrections to the energies of excited levels of light nuclei have 
been made in G4FermiBreakUp. 

After several fixes and updates have been introduced into G4StatMF, we 
have also achieved a rather good agreement with the results obtained with 
the FORTRAN SMM code. The validity of G4StatMF is confirmed in the 
stand-alone tests, which simulated the decay of hot nuclear systems with 
masses between Sn and Pb. We have also shown that the G4StatMF is able 
to describe central 129 Xe+ nat Sn collisions at 32 A MeV. The results of these 
tests and the proposed corrections were made available to the members of the 
Geant4 collaboration. The patches to G4StatMF were officially accepted and 
became a part of the Geant4 version 9.2 in December 2008. Now G4StatMF 
can be used to predict the production of various fragments and their isotope 
distributions in decay of hot medium-weight and heavy nuclei created in 
nuclear fragmentation reactions in extended media. 



As shown in our previous publications [15|, [16|, [17|, uM, the Light-Ion 
Binary cascade model is rather successful in describing the transport of light 
nuclei of therapeutic energies in water. The use of the Wilson abrasion 
model for such calculations is also validated by the present study. These 
models combined with the Fermi break-up model of Geant4 describe well the 
total dose and the yields of secondary fragments (H, Li, Be) from carbon 
beams. At the same time future work is needed to improve the description 
of He and B fragments by these models. We have to stress that not only 
theoretical models but also the measurements in extended targets should be 
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improved in the future, as currently available experimental data still have 
large uncertainties. 

The propagation of energetic Fe nuclei in water has been simulated using 
different models from the Geant4 toolkit. The initial stage of nucleus-nucleus 
collisions was described either by means of the Light-Ion Binary cascade 
model or by the Wilson abrasion model. As pointed out by developers of 
the cascade model, at least one of the colliding nuclei should be lighter than 
nitrogen. This condition is violated for Fe nuclei in water (H2O) that may 
be the reason of overestimated dose before the Bragg peak, possibly due 
to inaccurate estimation of sizes and excitation energies of Fe spectators 
by this model. There are no restrictions on the masses of colliding nuclei 
in the Wilson abrasion model, and this model seems more relevant to Fe 
fragmentation in water providing much better agreement with the depth- 
dose data before the Bragg peak. However, the dose behind the Bragg peak 
is underestimated by the abrasion model, possibly, due to underprediction of 
the yields of secondary protons and light nuclei. 

Our analysis shows that the Light-Ion Binary cascade model is not good 
enough to simulate fragmentation of heavy nuclei, despite the fact that its 
component, G4StatMF, is valid for stand-alone multifragmentation calcu- 
lations. We conclude that better models for the initial stage of nucleus- 
nucleus collisions are needed to describe the transport of heavy nuclei in ex- 



tended media by means of the Geant4 toolkit. As shown in Ref. [65], JQMD 
(JAERI version of Quantum Molecular Dynamics) model can be considered 
as a promising candidate. We believe that G4FermiBreakUp and G4StatMF, 
which have been validated in the present paper, can be coupled with such a 
model in the future. 
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